Overview

The mesoclim package provides functions to enable the mechanistic downscaling of climate data from coarse temporal and/or spatial resolutions to finer resolutions at which mesoclimatic processes, such as the effects of elevation, cold air drainage and coastal exposure, significantly affect how local weather conditions vary across the landscape. The package will permit the downscaling of climate data with a spatial resolution of 10s of kilometers to resolutions of under a single kilometer. It does not permit the modelling of microclimatic processes or below-canopy or below-ground conditions, but provides suitable data for subsequent microclimate modelling (refs to other packages).

The functions address key steps and applications of the downscaling work flow, namely:

{Something on using weather station inputs??}

Functions are also provided to carry out simple data checking, statistical summaries and graphing of climate datasets.

1 Downloading of coarse resolution global climate data.

A number of source-specific functions are provided for both the downscaling and processing of coarse-resolution climate data.Data sources include ERA5 Reanalysis data produced by the Copernicus Climate Change Service (ref), NCEP-NCAR and NCEP–DOE Atmospheric Model Intercomparison Project (Kanamitso et al 2002), and UKCP18 global and regional future climate estimates. {Others??}.

Example: downloading UKCP18 RCM data for the UK

UKCP18 data is available at different resolutions from various climate models that have been run to produce outputs for global, European and UK domains. They are the outputs of various comoinations of driving global climate models and regional or local

Data download requires a ceda username and password.

The function will download netcdf file(s) containing a decade of data that correspond to the model, collection, domain and time period requested. UKCP18 files contain a single variable and therefore multiple files will be downloaded

For further details see: {MetOff reference}

dir_ukcp18<-tempdir()
collection<-'land-rcm'
domain<-'uk'
member<-'01'
rcp<-'rcp85'
startdate<-as.POSIXlt('2018/01/01')
enddate<-as.POSIXlt('2018/12/31')

download_ukcp18(
    dir_ukcp18,
    cedausr,cedapwd,
    startdate,enddate,
    collection, domain, rcp, member,
    vars=c('clt','hurs','pr','prsn','psl','rls','rss','tasmax','tasmin','uas','vas')
)

list.files(dir_ukcp18)

Ancillary data requirements and sourcing

Digital terrain models

Additional data required for downscaling include digital terrain models (DTMs) at both the coarse resolution of the climate input data and at the fine resolution of the downscaling.

Where possible, the coarse resolution DTM should correspond to data that informed climate modelling. For example, this is made available for UKCP18 RCM data and can be derived from ERA5 using the geopotential variable. {An option in the UKCP18 and ERA5 download functions allows the downloading of these DTMs?}

Where the original DTM corresponding to input climate data is unavailable , if can be derived from various webservices and resampled to the correct projection and resolution.

The DTMs used in downscaling also act as a land/sea mask, and therefore sea cells should possess an NA value.

TO DO: * ADD download of UKCP18 rcm dem option to function and to ERA5 using geopotential!!!

Sea Surface Temperature

Sea surface temperatures are required for the corresponding time period and geography of the downscaling area. Functions are provided to download historic sea surface temperature from ERA5 and NCEP(?) data and future data estimated using UKCP18 regional models for the UK and NW Atlantic (add reference).

__ Something more on sourcing future sst data??__

# Download UKCP18 sea surface temp data

# ADD FUNCTION

Additional data may also be used to inform preprocessing and downscaling of particular climate inputs. For example, UKCP18 data commonly provides net rather than downward short and longwave radiation. For converting to the downward radiation required by downscaling, albedo values can either provided (from historic data sources) or constant proxy values can be used.

2 Preprocessing coarse resolution climate data

Defining area of interest

Unless we are wanting to analyse the whole area covered by downloaded climate data, it is best to define a broad area of interest to restrict data size and processing times. The area of interest is used to crop a coarse DTM of the same projection and resolution as the climate data to provide a template for data processing.

# Load UKCP18rcm DTM of UK
dir_data<-system.file('data-raw',package='mesoclim')
dtm<-terra::rast(file.path(dir_data,'orog_land-rcm_uk_12km_osgb.nc'))
# Crop to an area of interest
aoi<-terra::vect(terra::ext(-7.125,-2.875,49.375,51.625),crs='EPSG:4326')
aoi_e<-terra::project(aoi,terra::crs(dtm))
dtm<-terra::crop(dtm,aoi_e)

plot(dtm,main='DTM of climate data extent to be extracted')

Creation of standard inputs for downscaling

Different sources of climate and ancillary data can provide different variables using different SI units, in different formats and using different file naming conventions. Therefore, mesoclim provides several functions to convert and check source data into standard inputs for subsequent downscaling. Here we look at UKCP18 regional data as an example.

UKCP18 climate preprocessing

UKCP18 data preprocessing requires the conversion of net to downward short and longwave radiation. The conversion of shortwave radiation is calculated using an estimate of albedo at the same resolution as the climate data. Albedo data can either be provided or when not, constant land/sea albedo values are used.


# Preprocess UKCP18 data using constant albedo land / sea values
collection<-'land-rcm'
domain<-'uk'
member<-'01'
rcp<-'rcp85'
startdate<-as.POSIXlt('2018/01/01')
enddate<-as.POSIXlt('2018/12/31')

ukcp18rcm<-ukcp18toclimarray(dir_data, dtm,  startdate, enddate,
                             collection, domain, member)
#> Loading clt_rcp85_land-rcm_uk_12km_01_day_20101201-20201130.nc
#> Loading hurs_rcp85_land-rcm_uk_12km_01_day_20101201-20201130.nc
#> Loading huss_rcp85_land-rcm_uk_12km_01_day_20101201-20201130.nc
#> Loading pr_rcp85_land-rcm_uk_12km_01_day_20101201-20201130.nc
#> Loading psl_rcp85_land-rcm_uk_12km_01_day_20101201-20201130.nc
#> Loading rls_rcp85_land-rcm_uk_12km_01_day_20101201-20201130.nc
#> Loading rss_rcp85_land-rcm_uk_12km_01_day_20101201-20201130.nc
#> Loading tasmax_rcp85_land-rcm_uk_12km_01_day_20101201-20201130.nc
#> Loading tasmin_rcp85_land-rcm_uk_12km_01_day_20101201-20201130.nc
#> Loading uas_rcp85_land-rcm_uk_12km_01_day_20101201-20201130.nc
#> Loading vas_rcp85_land-rcm_uk_12km_01_day_20101201-20201130.nc
#> Using constant land and sea albedo values - assuming NA values in dtmc are sea!!

UKCP18 sea surface temperature preprocessing

ukcp18sst<-create_ukcpsst_data(dir_data,startdate,enddate,aoi,member)

plot(ukcp18sst, range=range(global(ukcp18sst,'range',na.rm=TRUE)))

Checking data inputs to downscaling

The resulting data structures of preprocessing can be checked to ensure there are no missing or unexpected values that may indicate a difference in the expected SI units or incomplete input datasets. This is particularly advisable if the inputs for spatial downscaling are not derived from one of the provided functions.

ukcp18rcm<-checkinputs(ukcp18rcm, tstep = "day")
#> [1] "Weather observations = 365"
#> [1] "Timesteps= 24 hrs, max= 24 hrs, min= 24 hrs"
#> [1] "Observations over 1 years, or 11.97 months, or 364 days."
#>              Min.    Mean    Max.
#> cloud       0.000  76.391 100.000
#> relhum     49.825  81.409  99.518
#> prec        0.000   2.904 134.241
#> pres       96.856 101.749 104.649
#> lwrad     186.724 317.669 395.398
#> swrad       3.090 139.858 352.741
#> tmax        2.823  14.267  29.030
#> tmin       -5.339  11.914  20.299
#> windspeed   0.003   5.036  22.359
#> winddir     0.001 168.808 360.000
#> elevation   8.017 114.508 284.199
#> [1] "Plotting spatial variation by month: red=max, green=mean, blue=min"

#> [1] "Plotting wind direction figures"

3 Spatial Downscaling

The projection, resolution and extent of the downscaled data is defined by providing a fine-resolution DTM that also acts as a fine resolution land/sea mask with sea cells having a NA value.

# Get fine and medium res DTMs from package
dir_datasets<-system.file('data',package='mesoclim')

load(file.path(dir_datasets,'dtmf.rda'))
load(file.path(dir_datasets,'dtmm.rda'))


dtmm<-unwrap(dtmm)
dtmf<-unwrap(dtmf)

Single step downscaling

Wrapper function allows downscaling in a single function call, with parameters defining which processes are captured by the downscaling.

Improvements

  • Include dtmf in spatialdownscale output.

  • Allow spatialdownscale output as spatrasters, packedspatrasters or arrays.

  • Include wind and temp height in spatialdownscale output

  • Write climdata to file function that converts to packed rast or arrays?

  • Read climdata file function?

# Interpolate SST data to required timesteps
sst<-spatial_interpNA(ukcp18sst)
sst<-time_interp(sst,ukcp18rcm$tme)

# Downscale one year of future climate data
t0<-now()
mesolizd<-spatialdownscale(ukcp18rcm, sst, dtmf, dtmm, basins = NA, cad = TRUE,
                           coastal = TRUE, refhgt = ukcp18rcm$tempheight_m, 
                           uhgt = ukcp18rcm$windheight_m, 
                           rhmin = 20, pksealevel = TRUE, patchsim = TRUE, 
                           terrainshade = FALSE, 
                           precipmethod = "Tps", fast = TRUE, noraincut = 0.01)
#> Downscaling wind...
#> Downscaling temperature...
#> Downscaling relative humidity
#> Downscaling pressure...
#> Downscaling SW radiation...
#> Downscaling LW radiation with terrain shading...
#> Downscaling precipitation...
#> Formatting output...
print(now()-t0)
#> Time difference of 6.321956 mins

# saveRDS(file = 'C:/Users/jm622/OneDrive - University of Exeter/Rprojects/mesoclim_local/Test_data/mesolizd.RDS', mesolizd)

# mesolizd<-readRDS('C:/Users/jm622/OneDrive - University of Exeter/Rprojects/mesoclim_local/Test_data/mesolizd.RDS')
# test<-lapply(mesolizd,terra::unwrap)

Display outputs

In tabular form for all times and locations:

smry_fun<-function(x) summary(as.vector(.is(x)))
rslt<-sapply(mesolizd[names(mesolizd)],smry_fun)
stats_df<-as.data.frame(t(round(rslt,3)))[,c('Min.','Mean','Max.')]
print(stats_df)
#>              Min.    Mean    Max.
#> tmin       -5.862  10.249  20.227
#> tmax        4.756  14.509  26.831
#> relhum     44.005  83.854 100.000
#> pres       94.888 101.003 104.588
#> swrad       6.218 132.752 351.011
#> lwrad     220.937 318.311 390.456
#> windspeed   0.034   3.950  24.487
#> winddir     0.000 173.102 360.000
#> prec        0.000   3.143  69.794

Or as selected rasters for days corresponding to spatial quantiles (eg days where spatial means are ata min, median and max values):

CHECK: Max diurnal temperature range plot shows grid effect.

for(var in names(mesolizd)){
  r<-terra::unwrap(mesolizd[[var]])
  names(r)<-rep(var,nlyr(r))
  plot_q_layers(r)
}


# Show spatial range in daily temperatures
plot(max(unwrap(mesolizd$tmax)-unwrap(mesolizd$tmin)),main='Max diurnal temperature range', font.main = 1, nc=1)

par(mar=c(1,1,1,1),cex.main=0.8, mgp=c(3,0.1,0))
layout(matrix(c(1,2,3,4,1,5,6,7,1,8,9,10),ncol=3),heights=c(1,3,3,3))

#par(mfrow=c(3,3))
plot.new()
text(0.5,0.5,"Spatial mean (green), max (red) and min (blue) by day of year",cex=1,font=1)

for(v in names(mesolizd)){
  r<-unwrap(mesolizd[[v]])
  plot_timestats_r(r,v,idx='doy',lgd=FALSE)
}

Multiple step downscaling

Spatial downscaling can be decomposed to the effects of specific physical processes, such as the effect of elevation, coastal or cold air drainage on temperature. Specific functions for downscaling particular effects can be applied without needing to undertake full downscaling using all processes.

Note: for individual functions, outputs may need further processing to ensure names are times are carried through.

Improvements

  • Allow option of spatrast or array outputs? Or enforce a standard spatrast output?

  • Allow passing standard data input list of climdata as well as individual variables

  • Allow variable winddir in windspeed downscaling - or warning message if input wind dir varies?

Pressure

pres<-presdownscale(ukcp18rcm$climarray$pres, dtmf, ukcp18rcm$dtmc, sealevel = TRUE)
names(pres)<-rep('pres',nlyr(pres))
terra::time(pres)<-ukcp18rcm$tme

plot_q_layers(pres)

Cold air drainage

Requires the estimation of drainage basins defined by topography which can be carried out separately from the estimation of cold-air drainage in downscaling.

# Basins - dtmf should have sea marked as NA
basins<-basindelin(dtmf,boundary=2)
plot(basins,main='basins')

For the chosen study area, cold air drainage has very little effect on local temperatures:

# Calculate using tmin to illustrate
ukcp18rcm$climarray$temp<-ukcp18rcm$climarray$tmin
tcad<-.tempcad(ukcp18rcm,dtmf,basins,refhgt=ukcp18rcm$tempheight_m)

names(tcad)<-rep('temp_cad',nlyr(tcad))
terra::time(tcad)<-ukcp18rcm$tme

plot_q_layers(tcad)

Temperature elevation downscaling

Variable lapse rates derived from temperature, humidity and pressure are applied to correct for elevation effects when downscaling. For daily temperatures this can be carried out for min and max daily values.

tminelev<-.tempelev(ukcp18rcm$climarray$tmin,dtmf,ukcp18rcm$dtmc,ukcp18rcm$climarray$relhum,ukcp18rcm$climarray$pres)
tmaxelev<-.tempelev(ukcp18rcm$climarray$tmax,dtmf,ukcp18rcm$dtmc,ukcp18rcm$climarray$relhum,ukcp18rcm$climarray$pres)

names(tminelev)<-rep('tmin_elev',nlyr(tminelev))
terra::time(tminelev)<-ukcp18rcm$tme

names(tmaxelev)<-rep('tmax_elev',nlyr(tmaxelev))
terra::time(tmaxelev)<-ukcp18rcm$tme

plot_q_layers(tminelev)

plot_q_layers(tmaxelev)


# Check dirunal range of elev downscaling
diurnaltmp<-tmaxelev-tminelev
names(diurnaltmp)<-rep('diurnal_t_range',nlyr(diurnaltmp))
plot_q_layers(diurnaltmp)

Wind

Windspeed downscaling aims to capture the effect of both elevation and the sheltering effects of topography that are in turn dependent on wind direction and height above ground.

windspeed<-winddownscale(ukcp18rcm$climarray$windspeed, ukcp18rcm$climarray$winddir, dtmf, dtmm, ukcp18rcm$dtmc, uz = ukcp18rcm$windheight_m)

names(windspeed)<-rep('windspeed',nlyr(windspeed))
terra::time(windspeed)<-ukcp18rcm$tme

plot_q_layers(windspeed)

Coastal effects

Uses sea surface temperatures, DTMs where sea is indicated by NA values and downscaled windspeed to estimate the effect of coastal effects. Coarse, medium and fine resolution DTMs are used to determine effects of upwind sea areas.

ADD: Guidance on selection of the extent of dtmm - how far should dtms extend from the downscaling area to adequately capture coastal effects??

sstf<-.resample(sst,dtmf)
tmincoast<-.tempcoastal(tminelev,sstf,windspeed,ukcp18rcm$climarray$winddir,dtmf,dtmm,ukcp18rcm$dtmc)
tmaxcoast<-.tempcoastal(tmaxelev,sstf,windspeed,ukcp18rcm$climarray$winddir,dtmf,dtmm,ukcp18rcm$dtmc)

names(tmincoast)<-rep('tmin_coastef',nlyr(tmincoast))
terra::time(tmincoast)<-ukcp18rcm$tme
names(tmaxcoast)<-rep('tmax_coastef',nlyr(tmaxcoast))
terra::time(tmaxcoast)<-ukcp18rcm$tme

plot_q_layers(tmincoast)

plot_q_layers(tmaxcoast)

All temperature effects downscaling

If only concerned with temperature downscaling, the wrapper function tempdownscale allows all relevant processes to be run as a single function.

ukcp18rcm$climarray$temp<-ukcp18rcm$climarray$tmin
tmin<-tempdownscale( ukcp18rcm,
                  sst,# can be coarse resolution
                  dtmf, # fine scale dem
                  dtmm , # medium re wider are dem
                  basins = basins, # basindelin() output or will calculate
                  u2 = windspeed, # windspeeds downscaled - or will calculate
                  cad = TRUE,
                  coastal = TRUE,
                  refhgt = ukcp18rcm$tempheight_m, # temp height
                  uhgt = ukcp18rcm$windheight_m # wind height
                  )

ukcp18rcm$climarray$temp<-ukcp18rcm$climarray$tmax
tmax<-tempdownscale(ukcp18rcm,
                  sst,# can be coarse
                  dtmf, # fine scale dem
                  dtmm , # medium re wider are dem
                  basins = basins, # basindelin() output
                  u2 = windspeed, # windspeeds downscaled
                  cad = TRUE,
                  coastal = TRUE,
                  refhgt = ukcp18rcm$tempheight_m, # temp height
                  uhgt = ukcp18rcm$windheight_m # wind height
                  )

plot_q_layers(tmin)

plot_q_layers(tmax)


diurnaltmp<-tmax-tmin
names(diurnaltmp)<-rep('Diurnal_temp_range',nlyr(diurnaltmp))
plot_q_layers(diurnaltmp)

Precipitation

The precipdownscale function provides various options and methods to downscale rainfall.

method<-'Elev'
precelev<-precipdownscale(
  ukcp18rcm$climarray$prec,
  dtmf,
  ukcp18rcm$dtmc,
  method = method,
  fast = TRUE,
  noraincut = 0.01,
  patchsim = FALSE,
  nsim = dim( ukcp18rcm$climarray$prec)[3]
)
names(precelev)<-rep('prec_elev',nlyr(precelev))
terra::time(precelev)<-ukcp18rcm$tme
plot_q_layers(precelev)


method<-'Tps'
prectps<-precipdownscale(
  ukcp18rcm$climarray$prec,
  dtmf,
  ukcp18rcm$dtmc,
  method = method,
  fast = TRUE,
  noraincut = 0.01,
  patchsim = FALSE,
  nsim = dim( ukcp18rcm$climarray$prec)[3]
)
names(prectps)<-rep('prec_tps',nlyr(prectps))
terra::time(prectps)<-ukcp18rcm$tme

plot_q_layers(prectps)


prec<-prectps

Longwave radiation downscale

Longwave downscaling accounts for the effects of terrain shading effect on skyview.

# simple downscaling
lwf<-.resample(.rast(ukcp18rcm$climarray$lwrad,ukcp18rcm$dtmc), dtmf, msk=TRUE)
# corrected for skyview  
svf<-.rta(.skyview(dtmf),dim(lwf)[3])
lwrad<-.rast(.is(lwf)*svf,dtmf)

names(lwrad)<-rep('lwrad',nlyr(lwrad))
terra::time(lwrad)<-ukcp18rcm$tme
plot_q_layers(lwrad)

Shortwave radiation downscale

Can be carried out without cloud patchiness or terrain shading…

swrad<-swdownscale(ukcp18rcm$climarray$swrad,as.POSIXlt(ukcp18rcm$tme, tz = "UTC"),dtmf,ukcp18rcm$dtmc,patchsim = FALSE,terrainshade = FALSE)

names(swrad)<-rep('swrad',nlyr(swrad))
terra::time(swrad)<-ukcp18rcm$tme
                       
plot_q_layers(swrad)

… or with cloud patchiness…

swradp<-swdownscale(ukcp18rcm$climarray$swrad,ukcp18rcm$tme,dtmf,ukcp18rcm$dtmc,patchsim = TRUE,nsim= length(ukcp18rcm$tme),terrainshade = FALSE)

names(swradp)<-rep('swrad_with_patchiness',nlyr(swradp))
terra::time(swradp)<-ukcp18rcm$tme
                       
plot_q_layers(swradp)